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SINGULAR PERTURBATION TECHNIQUES 
FOR REAL TIME AIRCRAFT TRAJECTORY 
OPTIMIZATION AND CONTROL 

♦ + 

Anthony J. Calise and Daniel D. Moerder 

Drexel University, Phiiadelphia, PA. 19104 

SUMMARY 


This study examines the usefulness of singular perturbation methods 
for developing real time computer algorithms to control and optimize 
aircraft flight trajectories. A minimum time intercept problem using F-8 
aerodynamic and propulsion data is used as a baseline. This provides a 
framework within which issues relating to problem formulation, solution 
methodology and real time implementation are examined. Theoretical 
questions relating to separability of dynamics are addressed. With respect 
to implementation, situations leading to numerical singularities are iden- 
tified, and procedures for dealing with them are outlined. Also, parti- 
cular attention is given to identifying quantities that can be precomputed 
and stored, thus greatly reducing the on-board computational load. Numeri- 
cal results are given to illustrate the minimum time algorithm, and the 
resulting flight paths. An estimate is given for execution time and 
storage requirements. 


Associate Professor, Dept, of Mechanical Engineering and Mechanics 
^Graduate Research Assistant 
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SECTION 1 


INTRODUCTION 


This report documents the derivation and evaluation of an on-line 
algorithm for minimum-time intercept. An approximation to the solution 
of the minimum-time intercept problem, resulting in a feedback control 
law, was derived using singular perturbation theory. The resulting 
control logic was evaluated using a model for the F-8 aircraft dynamics. 

The singular perturbation method is an order reduction procedure 
where a system’s dynamics are separated into fast and slow modes. There 
are two major benefits. First, higher-order problems can be approximated 
by a series of lower-order ones; and second, numerically "stiff” systems, 
having extreme differences in their modal behavior are solved on separate 
time scales. 

Applications of singular perturbation theory to flight mechanics 
problems have, to date, centered primarily on aircraft trajectory optimiza- 
tion. In the early seventies, a number of papers appeared in which energy 
state and reduced order modelling, with separate boundary layer correction 
terms, were employed to solve a number of flight control optimization 
problems [1-6]. Matched asymptotic expansions were employed somewhat later 
in solving the minimiam-time-to-climb problem, observing that the first 
order approximation matched well with the numerical solution obtained for 
the full order problem using a steepest descent technique [7]. Singular 
arcs have also been studied in the minimum-tirae-to-climb problem [8, 9]. 

In a more recent sequence of papers,, complete time scale separation 
has been applied in order to obtain feedback control solutions. Problems 
studied this way include the vertical plane minimum time and fuel problems 
and those involving weighted combinations of the two [10-15]. Feedback 



control of a missile in the horizontal plane was examined in [16]. The 
three-dimensional minimum;- time interception problem was examined by 
patching together two-dimensional subproblems in [17]. Under a NASA-sup- 
ported research effort that paralleled this effort, trajectory optimization 
in the long range, three-dimensional minimum-time intercept problem has 
been studied [18]. Singular perturbation theory has also been found 
useful in differential game problems [19-21]. The results derived 
in this study are mainly based on the approach in [12]. 

The research described in this report resulted in several extensions 
in the "state of the art" for the three-dimensional minimum-time interception 
problem; 

1. This work provides solutions for both long and short range cases, 
with short range cases being characterized by the lack of a 
maximum velocity cruise arc. 

2. A boundary layer matching criterion was derived such that the 
terminal boundary layer matching conditions could be calculated 
off-line, thus eliminating the problem of attempting to backwards 
integrate the terminal boundary solution on-line. This has been 
identified as one of the major stumbling blocks encountered in 
the application of singular perturbation theory [22]. Indeed, 
in the short range case, matching between initial and terminal 
layer arcs is the key issue to be resolved in solving the problem. 

3. The ordering and separation of energy and heading rate dynamics 
has been defined and justified for this problem, where range to 
interception is sufficiently large that optimal trajectories are 
dominantly characterized by climb and descent, rather than by 
turning. It should be noted that the ordering in this report is 
the reverse of that used in [1] and [10], and is due to the pre- 
sence of relative position dynamics in the model. 
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4. A detailed analysis has been made of issues pertaining to the 
separation of altitude and flight path angle dynamics. Optimiza- 
tion of altitude and flight path angle dynamics is a classical 
problem in flight mechanics, which to date had not yielded a 
solution suitable for on-line implementation. The essential 
problem is that these dynamics are highly coupled when considered 
apart from position and energy dynamics . It is shown that the 
dynamics, while not completely separable, can be approximated 

by singular perturbation methods with the inclusion of a penalty 
term on flight path angle in the performance index. 

5. The control algorithm was implemented as a flight director in the 
F-8 real time simulator at NASA Langley. Results of this evalua- 
tion are given in [23], and a NASA technical report is currently 
being prepared on this topic. The major issues relating to real 
time implementations have been addressed as part of the research 
effort. 

The organization of this report is as follows. Section 2 describes 
the problem formulation and summarizes the zero order singular perturbation 
solution. Section 3 gives a summary of the on-line control logic, describing 
the practical issues encountered during implementation of the control law. 
Section 4 describes the point-mass model of the F-8 aircraft employed in 
generating numerical results. Section 5 presents numerical results for this 
aircraft. Section 6 gives the conclusions and recommendations for future 
work. Detailed derivations of the control solution presented in Section 2 may 
be found in Appendices A and B. The first deals with position, energy and 
heading dynamics , and the second, with altitude and flight path angle 
dynamics. Appendix C presents the analysis relating to separation and 
ordering of energy and heading rate dynamics. Appendix D presents a means 
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by which the altitude and flight path angle dynamics can be analyzed on 
separate time scales, and gives a detailed analysis of the issues involved. 
Appendix E provides a description of and justification for a procedure 
for minimizing a Hamiltonian function with one unknown adjoint variable, 
employed several times in the analysis in Appendices A and B. A detailed 
description of the measures taken to suppress numerical singularities from 
the computed control solution is given in Appendix F, 
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SYMBOLS 


Drag Coefficient 

Zero Lift Drag Coefficient 

Lift Coefficient 

Slope of Lift Coefficient Curve, 1/rad 
Drag, N 

Drag for L = W, N 
Total Energy per Unit Weight, m 
Pseudo Cruise Energy Level, m 
Long Range Cruise Energy Level, m 

2 

Gravitational Constant, ra/s 

Hamiltonian 

Altitude, m 

Induced Drag Parameter 

Flight Path Angle Weighting Parameter in Cost Function 
Lift, N 

Vertical Lift Component, N 
Horizontal Lift Component, N 
Mach Number 
Mass, kg 

Energy State and Costate Eigenvalues, 1/s 
Heading State and Costate Eigenvalues, 1/s 

2 

Dynamic Pressure, N/m 
Range , m 

2 

Reference Area, m 
Thrust, N 
Time, s 
Time-To-Go, s 
Velocity, m/ s 
Weight, N 

Horizontal Position Variables, m 

Angle of Attack, rad 

Zero Lift Angle of Attack, rad 



Symbols 

6 

Y 

6 (.) 


£ 

n 

A 



(coat,) 

* Heading, rad 
= Flight Path Angle, rad 

» Perturbation Associated with a Particular Variable 
» Perturbation Parameter 
» Induced Drag Parameter 
- Line-of-Sight Angle, rad 
« Energy Costate, s/m 


X, “ Altitude Costate, s/m 

3 


X , X = Horizontal Position Costates, s/m 
^o 

X = Heading Costate, s/rad 


X = Flight Path Angle Costate, s/rad 

^4 



=» Bank Angle, rad 

3 

* Air Density, kg/m 
= Stretched Time Scale, s 
= Time Constants, s 


Subscripts 

c = Climb 

d = Descent 

D = Desired 

f = Final 

p = Relating to Proportional Vertical Lift 

max = Maximum 

min = Minimum 

T = Target 

o = Outer Solution Variable or Solution 

1, 2,3,4 = Boundary Layer Variables or Solution 

Superscripts 

(.) = Nominal Value 

o = Artificially Perturbed Variable 
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SECTION 2 

PROBLEM FORMULATION AND SOLUTION 

This section summarizes the problem formulation and the resulting 
singular perturbation solution. Derivations of the results summarized 
here may be found in Appendices A and B. Numerical results are given 
to illustrate the control solution for a variety of flight conditions. 


2.1 Problem Formulation 


The equations of motion are written in a horizontal, target centered 
coordinate frame: 


• 

X 

- 

V cos Y cos g 



(2.1) 

y 

m 

V cos Y sin 6 - cos y-j 



(2.2) 

£E 

» 

(T-D) V/W 



(2.3) 

2^ 
e e 

= 

L sin y/m V cos y 



(2.4)"^ 

z h 

= 

V sin Y 



(2.5) 

4 . 

e Y 

St 

(L cos y - W cos y)/^V 



(2.6) 

The variables 

1 in (2. 1-2. 6) are defined 

with the 

aid of Figure 1 

, where the 

subscript 

TTi^r 

' is used to designate the 

target. 

These equations 

are valid 


for a flat earth, with thrust (T) directed along the flight path and 
constant weight. Drag (D) is assumed to have a conventional parabolic form 

D = qs Cjj = qs(Cp + n CI^) (2.7) 

a 

which can also be written as 

D = qs (Cjj + KL^/qs), q = pv^/2 (2.8) 

where q is the dynamic pressure, p is the air density and 

K = T)/C^ (2.9) 

a 

L =* qs » qs(Cj^ a) (2.10) 

a 

The variable E is the total aircraft energy (kinetic and potential) per unit 
weight . 


+ 


In this report. 


w/xyz is to be interpreted as w/(xyz). 
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E 


h + V^/2g 


( 2 . 11 ) 


where h is the aircraft altitude. The control variables are aircraft lift 
(L) , bank angle (y) and thrust (T) . The objective is to find the controls 
L, y, T that minimize 

^f 2 

J * / (1 + k sin y)dt , k ^ 0 (2.12) 

o 

where k » 0 for minimum time. The minimization is subject to the following 
state and control variable constraints: 


L < W G 
— max 


L < qs C_ a 
— ^ L max 
a 


T . (h, V) < T < T (h, V) 

min ' — “ max ’ 


q < q , M < M (h) 
— ^max — max 


(2.13) 

(2.14) 

(2.15) 

(2.16) 


where G is the maximum load factor, ot is the stall angle of attack, 
max max 

T . and T are the minimum and maximum thrust levels that are functions 
min max 

of aircraft altitude and velocity (V) . The boundary conditions are such that 
the initial aircraft state is fully specified, and we require 

x(t^) =» y(t^) » 0 , h(t^) = (2.17) 

for intercept, when h^(t^) is taken as the projected target motion in 
altitude 

hi(tf) = h^(0) + (V^ sin Y^) t^ (2.18) 

The objective here in using singular perturbation methods is to. approxi- 
mate the open loop optimal control solution with a near-optimal control 
solution in feedback form. Towards this end, the equations of motion in 
(2. 1-2. 6) have been scaled by powers of e, which imply that a natural separa- 
tion in the system dynamics exists. Ideally, one would like to identify e 
with small physical system parameters. This can be done with varying degrees 
of success by writ ting the equations of motion in a non-dimensional form. 

An example is given in [15]. The ordering selected in (2. 1-2. 6) is based 
on an understanding of aircraft dynamics, experience with problems in 
trajectory optimization, and the earlier results of researchers in this 
field. The approach here is to seek a solution for £“1 by an expansion 
about e»0. VThile this departs from the spirit of asymptotic expansions. 
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where e is regarded as "sufficiently small", the accuracy of the resulting 
solution depends more on the degree to which the dynamics are separated. 

The particular ordering selected here can be argued on a physical 
basis. Long range optimal trajectories are generally made up of climb, 
cruise and descent arcs. The climb and descent arcs can be considered as 
boundary layers needed to satisfy initial and terminal constraints on 
energy and altitude, during which energy is first increased (to the cruise 
energy) and then decreased (to satisfy the terminal altitude constraint). 
Initial transitions to the optimal climb path and the optimal heading for 
intercept take place on a much shorter time scale, in comparison to the 
time needed to gain or lose energy. Selecting 6 dynamics as slower than 
h and y d 3 mamics allows for high and low speed yo-yo maneuvers during the 
initial turn at large heading errors. This is illustrated in the results 
of Section 2.4. A detailed analysis of the ordering of E, B and h, y dy- 
namics is given in Appendices C and D- 

2.2 Outer Solution 

In the outer solution, the controlled aircraft is assumed to be trav- 
elling on a fixed course at a constant speed, as can be seen by letting 
e 0 in (2. 1-2. 6). The problem is reduced to optimal intercept in the hori- 
zontal plane. The state variables are x and y, and the controls are 3, 
h and E. In order to satisfy the intercept requirement (see Fig. 1), we 
must have 


Vsin(B-X) ® cos cos X (2.19) 

or, in other words, there is no relative motion allowed perpendicular to the 
horizontal plane line-of-sight axis. The optimal controls h^ and E^ are 
determined as 


h , E = arg max (V) (2.20) 

° ° h,E 

where the maximization takes place subject to the constraints in (2.13-2.16) 
and subject to 


T 

o 


M = 0, Y 


0. L 


W 


where D is drag for L = W 
o 

qs + KW^ /qs 
o 


and 

q - P(h,)V^/2 


V 

o 


■'28«o-"o> 


( 2 . 21 ) 

( 2 . 22 ) 

(2.23) 
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These constraints arise from setting e->0 in (2. 2-2,6). The cruise point 
(2.20) for the F-8 aircraft is displayed in Figure 2, superimposed on the 
aircraft’s flight envelope. It should be noted that, in a higher perfor- 
mance aircraft for which the q or Mach boundaries would be encomtered 
while T D, the solution would lie at the intersection of the constraint 
boundary with the zero energy rate boundary. This is illustrated for an 
F-4 aircraft in [25]. 

The optimal cruise heading (3^)is computed using (2.19): 

- sin”^ {V^ cos cos X/V} + X (2,24) 


The costates and , associated with the horizontal position dynamics 
in the outer soSution, while not explicitly appearing in the outer control 
solution, are used in subsequent boundary layer solutions. These take the 
form 


= -cos6^/(V^-V^ cos cos 3^) (2.25) 

o 

X = -sin3^/(V^-V^ cos cos 3^) (2.26) 


It should be noted that the cruise solution for h and E is indepen- 
dent of target motion and intercept geometry. This allows these quantities 
to be calculated off-line and stored. The only outer solution calculations 
performed on-line are (2.24-2.26). 

2.3 First Boundary Layer Solution 

The first boundary layer solution addresses energy dynamics. The 
constraints 


Pi - 0, = 0, - W (2.27) 

in addition to (2.13-2.16), arise when the time transformation x = t/e 
is introduced and we let e -►O. The controls are T, h, and 3. The optimal 
heading (3^) is identical to that for the outer solution. Since T appears 
linearly in the dynamics, we have 



, X < 0 

(2.28) 

■«ln 

. X > 0 

^1 

(2.29) 
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where is the energy costate 
essentially correspond to climb 
yields 


(T -D )V 
^ r TOax O T 
arg mln{ } 


variable* The solutions in (2«28) and (2*29) 
and descent. Optimization with respect to h 


E - E 

current 
T > D 


(2.30) 


for ascent, and 

-(T . -D )V 

\ = arg mn{ — } 

h o 

for descent . 

The climb path to cruise for the F-8 is superimposed on the aircraft 
flight envelope as a bold line in Figure 2* The optimal descent path is 
along the boxmdary. The expression for the first boundary layer 

costate is 


(2.31) 


E = E 
T < 


current 

D 




(2.32) 


where H (E, h,) is the outer solution Hamiltonian evaluated at the first 
o 1 

boundary conditions: 


H (E,h, ) = {X V cos B + X (VsinB - V cosy™) + 1} 
o 1 X y i i 

o o 

Since the solution for h^(E) is independent of target motion, it can 
be precomputed and stored as a function of E. Only the costate variable 
in (2.32) is computed on-line. 

2.4 Second Boundary Layer 

This boundary layer is obtained by introducing the time transformation 
2 

T =» t/e and letting e-K>. This results in the constraints; 

= 0 , =■ 1^2 + (2.34) 


(2.33) 

E = E 

current 
h = h. 


where L is the total lift and 1^2 the horizontal lift component. The 
controls in this boundary layer are T, h, and ^ 22 * Assuming that all turning 
takes place near the initial time where X <0, the optimal thrust is 


T- = T (h«, V.) 
2 max 2* 2' 


(2.35) 
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Figure 2. Flight envelope for the F-8 aircraft. 


( 2 . 36 ) 


where is defined by 


h 


2 


erg min {-p/KV IL(E,h,3)} 
h ^ 


E 

6 


E 

current 

6 

current 


H^(E, h, B) is the first boundary layer Hamiltonian evaluated at current 
values of E, h and 3. It is expressed as 

H^(E, h, “ ^ 2 ^ V cos B + A [Vsin3 - cos y^] + 

^o ^o 

X (T-D^)V/W + 1 (2.37) 

El o 

The solution for L ^2 is analytic, and is given by 

L ^2 * /-“qs W H^(E,h,3) * sign (B^-B) (2.38) 

After performing the minimization in (2.36), the heading costate variable 
is computed using 

- -2 H-(E,h,3) mV/L..j (2.39) 

®2 ^ h - h2 

The calculations in (2.36-2.39) must be performed on line. To accelerate 
the minimization in (2.36), which is performed each time the control 
solution is updated, the solution from the previous time instant is used 
as a starting point. 

Numerical results for the F-8 aircraft are given in Figures 3 and 4. 

These display L^o and h« , respectively, as functions of A3 = B -B for several 
zz z o 

values of E. It can be seen from the figures that, at all energy levels, 

L 22 Q ^2 ^ ^1 A3“^. This type of asymptotic behavior is necessary 

for a valid singular perturbation solution. Note also that at higher energies, 
h 2 does not digress very far from and that 1^22 even for large heading 
errors, is significantly below This indicates that, at these higher 

energies, the zero-order solution attempts to preserve the combination of 
energy rate and closure rate from the first boundary layer solution. At 
low energy levels the situation is considerably different, with emphasis 
placed on rapidly reducing heading error. Horizontal lift is saturated before 
A3 reaches 1.2 rad. For larger heading errors, h^ increases so that the 
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4696 m 



Figure 3. Normalized second boundary layer 

horizontal lift as a function of heading 
error for several energy levels. 







aircraft is driven to the comer velocity for that energy level. At a 
midrange energy, priorities are more evenly mixed. The solution for h 2 
remains close to h^ longer than in the low energy case. Also note that 
L 22 take on both saturated and nonsaturated values after encountering 

the maximum lift constraint boundary* These results illustrate the trade-off 
between energy rate, closure rate and turn rate that takes place in the 
second boundary layer solution. 


2.5 Third Boundary Layer Solution 

The third boundary layer addresses the altitude dynamics. It arises 

3 

from introducing the time transformation t *» t/e and letting e 0. The 
vertical lift this boundary layer is constrained to be 

^13 * ^ ^ (2-40) 

2 

Due to non-zero values of y, the penalty term k sin y in (2-12) has an 
effect on the solution for the third and fourth boundary layers. The 
controls here are horizontal lift (^ 23 ^ zero-order solution 

for L 22 is given by 

L 23 = min {L 2 ^ 3 ^, L 22 /COS Y) (2.41) 


where 


^2 max 



max 



(2.42) 


It can be seen that L 2 ^ L 22 as y 0 which is the constrained value 

for y in the second boundary layer. The expression for y^ is 

y^ = arg max {sin y/H 2 (h,E ,6 ,y) } •sign (h 2 "h) (2-43) 

y 


where H 2 i® second boundary layer Hamiltonian evaluated at the current 

conditions for its arguments: 

H 2 (h,E, 3 ,y) “ (X^ cos3 + X sin3)V cos Y - X cos y^ 

o ^o ^o 

+ (T-D)V/W + g/WV cos Y + 1 + ksinS (2.44) 

^1 ^2 


Numerical results for y^ as a function of altitude error for zero 
heading error is shown in Figure 5 . The effect of the parameter k is 
also shown. Note that y^ 0 as h 2 h^ and that, as might be expected, 

increasing k decreases the magnitude of y^ at all flight conditions. 

The costate for this boundary layer X, is determined from 
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(2.45) 


- - H 2 (h,E,e,Y 2 )/V sin 


2,6 Fourth Boundary Layer Solution 

In the fourth boundary layer, the vertical and horizontal components 

of lift (L^^ and respectively) are refined to reflect the flight path 

angle dynamics. As long as L $ L , the lift components are: 

” max 

- W cos Y - H 3 (h,E, 6 ,Y)qsW/X^ KV^ cos y (2.46) 


1*24 * gqs/2Xg^ KV^ cos y (2.47) 

These are used to define the final lift and bank angle commands: 

L - *^^4 + L 24 (2.48) 

y = tan ^ (L^^/L^^) (2.49) 


If L in (2 . 48) exceeds 
tan y = (X /X 

In this formulation X 

Y 


L - we set L = L and obtain an expression for y 
max max 

cos y) • sign (6^-0) (2.50) 

, the costate for y, is evaluated as a root of 
4 


+ BX + C =■ 0 

(2.51-a) 

^4 

g^(L - cos^y)/V^ 

max 

(2.51-b) 

2^ cos Yg/V 

(2.51-c) 

y)^ - 

P 2 max 

(2.51-d) 

H,-X„ (L^ -W^cos^Y)KV/qsW + X, V sin y 

1 E- max n 

(2.51-e) 


It is demonstrated in Appendix B that (2.51-a) will always have real roots 

of opposite sign; thus X is chosen such that 

■^4 


sign {X } = - sign (Yo-y) 
Y4 


(2.52) 
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Given this value for X , u is computed using (2.53) placing ]x in the 
quadrant appropriate to sign {y 3 “Y}« 

It should be noted that the arbitrary separation of h and y dynamics 
in the third and fourth boundary layers fails to account for the coupling 
that naturally exists between these states. A method is given in Appendix 
D for choosing k in (2.12) so that this problem is alleviated. 
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SECTION 3 


FEEDBACK IMPLEMENTATION 

This section describes the feedback implementation of the control 
solution formed from the singular perturbation outer and boundary layer 
solutions described in Section 2. Five topics are covered: organization 

of climb and descent legs, an alternate proportional vertical lift scheme, 
thrust and lift control during descent, avoidance of singularities in 
the control solution, and the overall organization of the actual feedback 
implementation • 

3.1 Climb and Descent Legs 

Long range intercept trajectories, ignoring initial turning and 

other transients, have three stages. The first stage is a climb to 

cruise at the long range optimal cruise energy (E*) - The second stage 

is a cruise leg. The third stage can take one of two forms. If the 

target altitude is below the altitude for long range cruise, it is a 

descent leg. It is important to note that climb and descent in this 

report refer to gain and loss of energy - not altitude . For example, 

altitude decreases during a portion of the climb profile. If the 

target altitude is above the long range cruise altitude, the terminal 

stage is a zoom climb (constant energy altitude gain) maneuver. 

Short range intercepts are defined as occurring when the intercept 

range is less than the range required to fly a long range climb and descent. 

In this case, the optimal trajectory would consist of climb and descent 

* 

that meet at an energy level less than at an altitude and velocity on 

the zero-energy-rate boundary for T * T^^^, or on the dynamic pressure or 

Mach constraint boundary (see Fig. 1). Henceforth, these lower energy 

points will be referred to as pseudo cruise points, E^. 

An important element in the control design is the decision logic 

required for determining whether an intercept path is long or short range. 

In the short range case, the logic must select a pseudo cruise energy 

level such that the horizontal range for climb and descent matches the 

* 

predicted intercept range. In the long range case, E is used and a 

o 

cruise leg is inserted to march the predicted range to intercept. Descent 
is initiated when the horizontal range for descent from E^ matches the 
predicted range to intercept. 
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As was mentioned in Section 2 , the climb, cruise and descent solutions 

are independent of intercept geometry and target parameters. The calctila- 

tions are performed off-line, and the results stored in the form h (E,E ) 

c o 

for climb and for descent, where E^ represents selected pseudo cruise 

energy levels. For the long range case we have h^(E,E^). Several climb 
altitude profiles for the F-8 aircraft are illustrated in Figure 6. The 
descent profile is independent of the selected cruise energy because 
satisfaction of (2.31) for h^ results in values on the dynamic pressure 
boundary. 

The time, t^(E,E^) , and horizontal distance, r^(E,E^), required to 
climb from E to E^ were determined by computing the integrals: 

E 

t^(E,E^) - /g° (l/E)dE (3.1) 

E 

r^(E.E^) - /g° (V^/E)dE (3.2) 


where 


-,/(E-hj^(E,E^))2g 


(3.3) 


and E is the energy rate computed at h^(E,E^) . Tabular data for h^(E,E^), 
t^(E,E^) and r^(E,E^) is presented in Table 1 for E^ and several pseudo cruise 
energies. 

The expressions used in calculating altitude and range for descent 
are slightly different from those used in calculating altitude and range 
for climb, since in descent, flight path angle is too large to be ignored. 

We have: 


t (E) = (l/E)dE 
o 

r^(E) = /\ (r^/E)dE 

where 


= /[E-h^(E)]2g 
= sin"^{(dh^/dE)E/Vj_} 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 
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TABLE 1 


F-8 CLIMB ALTITUDE, TIME AND RANGE 
AS FUNCTIONS OF ENERGY 


ASCENT VARIABLES 

« CRUISE 

ENERGY- 4696 

• 05 H 

ENER6Y(N) 

HCMl 

TIME(SEC) 

RAN6ECR} 

279.45 

0.0 

55.0058 

10342.0 

1383.60 

0.0 

34.4565 

7888.8 

2487.75 

0.0 

21.4960 

5389.9 

3591.90 

30.48 

10.4930 

2720.8 

4696.05 

820.36 

0.0 

0.0 

5800.20 

0.0 

0.0 

0.0 

6904.35 

0.0 

0.0 

0.0 

8008.50 

0.0 

0.0 

0.0 

9112.65 

0.0 

0.0 

0.0 

10216.80 

0.0 

0.0 

0.0 

11320.95 

0.0 

0.0 

0.0 

12425.10 

0.0 

0.0 

0.0 

13529.25 

0.0 

0.0 

0.0 

14633.39 

0.0 

0.0 

O.Q 

15737.55 

0.0 

0.0 

0.0 

16841.70 

0.0 

0.0 

0.0 

17945.85 

0.0 

0.0 

0.0 

19049.99 

0.0 

0.0 

0.0 


ASCENT VARlABLESt CRUISE 

ENERGY= 6904 

.35 M 

ENERGYCM) 

HCM) 

TIME(SEC) 

RANGE(R) 

279.45 

0.0 

76.8073 

16350.4 

1383.60 

0.0 

56.2579 

13697.2 

2487.75 

0.0 

43.2975 

11398.3 

3591.90 

30.48 

32.2945 

8729.2 

4696.05 

822.96 

21.8001 

6009.3 

5800.20 

1615.44 

11.1350 

3118.5 

6904.35 

2477.75 

0.6 

0.0 

8008.50 

0.0 

0.0 

0.0 

9112.65 

0.0 

0.0 

0.0 

10216.80 

0.0 

0.0 

O.Q 

11320.95 

0.0 

0.0 

0.0 

12425.10 

0.0 

0.0 

0.0 

13529.25 

0.0 

0.0 

0.0 

14633.39 

0.0 

0.0 

0.0 

15737.55 

0.0 

0.0 

0.0 

16841.70 

0.0 

Q.O 

0.0 

17945.85 

0.0 

0.0 

0.0 

19049.99 

0.0 

0.0 

0.0 
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TABLE 1 (CONTINUED) 


ASCENT VARIABLES! CRUISE 

ENERGY* 9112* 

65 H 

ENER6Y(M> 

HCH> 

TIHE<SEC) 

RANGECM) 

279,45 

0*0 

107*9332 

25825*6 

1383*60 

0*0 

87.3839 

23372*3 

2487.75 

0*0 

74*4234 

20873.4 

3591*90 

30*48 

63*4204 

18204*4 

4696*05 

822*96 

52*9260 

15484*5 

5800*20 

1615*44 

42*2609 

12593*7 

6904*35 

2407*92 

31*1056 

9437*3 

8008*50 

3139*44 

17*7964 

5471*5 

9112*65 

4065*64 

0.0 

0*0 

10216*80 

0*0 

0*0 

0*0 

11320*95 

0*0 

0*0 

o.d 

12425*1 0 

0*0 

0*0 

0.0 

13529*25 

0*0 

0*0 

0.0 

14633*39 

0.0 

0*0 

0*0 

15737.55 

0*0 

0*0 

0*0 

16841*70 

0*0 

0*0 

0*0 

17945*85 

0*0 

0*0 

0*0 

19049*99 

0*0 

0*0 

0*0 


ASCENT VARIABLES, CRUISE 

ENERGY= 11320* 

95 M 

ENERGY(M) 

H(M) 

TIME(SEC) 

RANGE(H> 

279.45 

0.0 

194*4020 

54474.5 

1383*60 

0*0 

173.8526 

52021.2 

2487*75 

0*0 

160.8922 

49522.3 

3591*90 

30*48 

149.8892 

46853.3 

4696*05 

822*96 

139.3948 

44133.4 

5800*20 

1615*44 

128.7297 

41242.6 

6904*35 

2407.92 

117.5744 

38086.2 

8008*50 

3139.44 

104.2652 

34120.5 

9112*65 

3870*96 

84.9084 

28070.4 

10216*80 

4541 *52 

51.5031 

17163.0 

11320*95 

5626*88 

0.0 

0.0 

12425*10 

0*0 

0.0 

0.0 

13529*25 

0.0 

0.0 

0.0 

14633*39 

0*0 

0.0 

0.0 

15737.55 

0.0 

0.0 

0.0 

16841*70 

0*0 

0.0 

0.0 

17945*85 

0*0 

0.0 

0.0 

19049*99 

0*0 

0.0 

0.0 
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TABLE 1 (CONTINUED) 


ASCENT VARIABLES* CRUISE ENERGYs 13529.25 M 


ENERGYCN) 

H(H) 

279.45 

0.0 

1383. GO 

0.0 

2487.75 

0.0 

3591.90 

30.48 

4696.05 

822.96 

5800.20 

1615.44 

6904.35 

2407.92 

8008.50 

3230.88 

9112.65 

4358.64 

10216.80 

5486.40 

11320.95 

6096.00 

12425.10 

7193.28 

13529.25 

7034.38 

14633.39 

0.0 

15737.55 

0.0 

16841.70 

0.0 

17945.85 

0.0 

19049.99 

0.0 


ASCENT VARIABLES* CRUISE 

ENERGY(M) 

H<M> 

279.45 

0.0 

1383.60 

0.0 

2487.75 

0.0 

3591.90 

30.48 

4696.05 

822.96 

5800.20 

1615.44 

6904.35 

2407.92 

8008.50 

3566.16 

9112.65 

4663.44 

10216.80 

5791.20 

11320.95 

6888.48 

12425.10 

7985.76 

13529.25 

8717.28 

14633.39 

9784.08 

15737.55 

8370.48 

16841.70 

0.0 

17945.85 

0.0 

19049.99 

0.0 


TIHE(SEC) 

RANGE(H) 

237.1549 

67539.1 

216.6056 

65085.9 

203.6451 

62586.9 

192.6421 

59917.9 

182.1478 

57198.0 

171.4827 

54307.2 

160.3274 

51150.8 

147.3731 

47332.3 

131.9121 

42741.9 

114.7642 

37635.2 

92.6095 

30741.9 

62.4630 

21150.7 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


ENERGYS 15737 

.55 H 

TIME(SEC) 

RAN6E<H> 

275.9463 

78177.2 

255.3971 

75724.0 

242.4366 

73225.1 

231.4336 

70556.0 

220.9393 

67836.1 

210.2742 

64945.3 

199.1189 

61789.0 

187.1362 

58435.7 

173.9817 

54709.5 

159.2435 

50510.3 

142.5676 

45719.5 

123.3140 

40147.4 

99.3018 

32953.4 

68.3067 

23472.8 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 
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TABLE 1 (CONCLUDED) 


ASCENT VARIABLES* CRUISE 


ENERGYCM) 


279.45 

0.0 

1383.60 

0.0 

2467.75 

0.0 

3591.90 

30.48 

4696.05 

822.96 

5800.20 

1615.44 

6904.35 

2407.92 

8008.50 

3596.64 

9112.65 

4846.32 

10216.80 

5974.08 

11320.95 

7071.36 

12425.10 

6168.64 

13529.25 

9265.92 

14633.39 

10027.92 

157-37.55 

11125.20 

16841.70 

10271.76 

17945.85 

10668.00 

19049.99 

10636.09 


ENERGY= 190^»9.9S H 


T.IHE(SEC) 

RANGE(H) 

628.5549 

210690.6 

608.0056 

208237.4 

595.0452 

205738.5 

584.0422 

203069.4 

573.5481 

200349.5 

562.8831 

197458.7 

551.7278 

194302.4 

539.8125 

190987.3 

527.1221 

187502.6 

513.2688 

183662.9 

497.5857 

179272.1 

479.4856 

174160.9 

458.2539 

168122.7 

431.6506 

160315.6 

397.0964 

149984.7 

343.9612 

132477.7 

244.0769 

95683.5 

0.0 

0.0 


where E is the energy at initiation of descent. Tabular values for h,(E), 

d 

t^(E) and ^^(E) are provided in Table 2. Climb times and ranges from 
the current energy to the cruise energy are obtained by interpolating 
and differencing the values in Table 1. A similar procedure is used 


for descent using Table 2. For example, 

td(E^»V " - "d<V ^3.9) 

rd(Ej,,Ef) - r^(Ef) - r^(E^) (3.10) 

In general, E^ is not known a priori and must be determined such that 
h(t^) =* h^(t^). A terminal constraint must be satisifed: 

hd(Ef) = sin (3.11) 

where t^^ is the estimated time remaining until intercept. Referring to 
Figure 1, we initiate descent when 

rd(Eo*Ef) ^ R + (V^ cos sin X)t^(E^,E^) (3.12) 

is satisfied, where R is the current horizontal range. 

3.2 Proportional Vertical Lift 


An option was included in the control logic for stopping the 
singular perturbation solution after the second boundary layer, and using 
a suboptimal proportional control for vertical lift. We first define a 
desired flight path angle 

= (h2-h)/T^V +E(dh^/dE)/V^ , - /TE-h^)2g (3.13) 

The second term in (3.13) is an approximation to the flight path angle for 
following the first boundary layer climb path. The proportional vertical 
lift (L^p) is computed based on a desired flight path angle rate proportional 
to (Yjj-y) : 

Yjj = (Yj)-Y)/t2 = (Lj^p-WcosY)/mV (3.14) 

Solving for we have 

Lip = mV(Yjj-Y)/t 2 + Wcosy (3.15) 
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TABLE 2 


F-8 DESCENT ALTITUDE, TIME AND RANGE 
AS FUNCTIONS OF ENERGY 


DESCENT VARIABLES* CRUISE ENERGY^ 19100.07 H 


ENERGYCM) 

K(N) 

3586.93 

0.0 

4499. A7 

677.40 

5412.00 

1331.60 

6324.54 

1976.64 

7237.08 

2604.15 

8149.62 

3212.69 

9062.16 

3814.46 

9974.69 

4390.52 

10887.23 

4957.39 

11799.77 

5506.95 

12712.30 

6030.09 

13624.85 

6553.26 

14537.38 

7055.07 

15449.92 

7531 .21 

16362.46 

8006.49 

17275.00 

8463.37 

18187.54 

8896.07 

19100.07 

9317.89 


TIHECSEC) 

RANGE(H) 

357.6738 

107306.4 

270.8616 

83918.9 

198.3982 

63759.6 

143.5021 

47992.5 

108.2608 

37554.1 

87.5223 

31231.0 

74.0457 

27015.2 

63.7786 

23719.3 

55.3226 

20935.3 

48.1671 

18520.0 

41.5890 

16232.8 

35.1294 

13919.1 

28.8112 

11585.8 

22.6819 

9252.1 

16.7482 

6929.9 

10.9957 

4614.5 

5.4142 

2303.1 

0.0 

0.0 




3-3 Thrust and Lift Control During Descent 

In the ideal case of a fully optimal control solution, there would 
be insignificant maneuvering and throttle variation during descent. There 
is, however, significant turning in the zero order solution implemented. 
This is primarily due to two factors. First, since the aircraft follows 
the dynamic pressure constraint boundary during most of the descent, and 
since the flight-path angle is non-zero, the intercept heading changes 
from the optimal cruise heading value. Because of this, it is necessary 
to update the intercept heading during descent using (2.24) and the hori- 
zontal component of aircraft velocity. Second, target maneuvers that 
occur after the initiation of descent necessitate heading changes. The 
former problem could be greatly reduced by correcting the outer solution 
to first order in E, in a manner similar to the procedure followed in [16 ] . 

In order to insure intercept under all conditions, it is necessary 
to modulate both thrust and descent. In the case ^^^22 

maximum lift should always be used during descent in maintaining the 

intercept heading, due to the fact that in (2.32) becomes positive 

^1 

during descent. In practice, a proportional control law is used 
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for small heading errors. 

Thrust modulations are used to control rate of descent to ensure that 
h«h^ at intercept. Two correction terms are introduced: 

T = T^ + 6 T^ + 6 T^ (3.21) 

where T, is the nominal descent thrust from (2.28, 2.29). The second 
a 

term corrects for the fact that L does not equal W during descent, since 
this is assumed in generating the t^ and r^ tabular data. Thus, 6 T^ com- 
pensates for the increased drag due to lift variations. 

6 T^ =• K(L^-W^)/qs (3.22) 

The second component compensates for the current mismatch ( 6 R)in range, 
where from (3.12), replacing by E 

5R » R + (V^ cos sinX)t^(E , E^)-r^(E ,E^)cos (B^-\) (3.23) 

A proportional control law was derived, defining 

6 R = -K^ 5R (3.24) 

Noting that 

6 R = (dr ,(E,E.) /dE)cos(3 -X) 6 E (3.25) 

at o 

and 

6 E = 6 T 2 V/W (3.26) 

one can solve for 6 T 2 as 

6 T 2 = K^W5R/Vcos(3Q-X)(dr^ (E,E^)/dE) (3.27) 

In order to allow for thrust variation, T . was set equal to T .,/2, 

* min ^ mil * 

where ^ is the military thrust level. 

It should be noted that a portion of the descent path calls for 

T^ = ^max* Referring to Figure 8 , note that upon initiating descent, 

for E > E , the commanded altitude on the descent path is such that 
os* ^ 

V, (E ) > V (E ) , As shown in Appendix A, this means that X_ remains 
1' o o o ’ E, 

negative until V^(E) < ^q(Rq)> which from (2.28) implies that T =» 
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3.4 Summary of Control Calculations 

A summary of the control calculations including proportional vertical 
lift is presented in Figure 9, The organization of the additional calcu- 
lations for the full zero-order solution is given by Figure 10. All 
on-line calculations are referenced by equation numbers. The purpose of 
the range matching block is to establish the proper cruise or pseudo cruise 
energy level, E^. During climb, the costates are calculated on-line and 
h^(E,E^) is taken from tabular data. All turn parameters are calculated 
on-line. During descent, h^(E,E^) is drawn from tabular data, and thrust 
and lift are calculated as described in Section 3.3. 

3.5 Avoidance of Numerical Singularities 

Numerical difficulties evidenced by discontinuities in the control 
solution were encountered when the aircraft altitude, heading and/or 
flight path angle approached their optimal values for the second and 
third boundary layer solutions. These discontinuities occurred when 
certain functions approached an indeterminate (zero over zero) form as 
the optimum state values were approached. For example, the argtiment 
being minimized in (2.36) approaches an indeterminate form as heading 
error approaches zero. The corrective measures took the form of first-order 
Taylor series expansions and approximations taking advantage of the 
asymptotic character or the boundary layer solutions. A more detailed 
description is given in Appendix F. 
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Figure 9. Summary of control calculations with 
proportional vertical lift. 
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Figure 10. Zero order lift and bank angle 
calculations. 




SECTION 4 


POINT-MASS MODEL OF THE F-8 AIRCRAFT 

This section describes the procedure by which the F-8 was modelled 
at trim conditions and presents the results of that effort. The inter- 
cept control law has been built on the assumption that the aircraft 
behaves as a point mass. The aerodynamic model takes the form: 

D * qs Cjj + ^ S (^*1) 

0 

L - qs C^ (4-2) 

a 

Thus, Implementation of this model for the F-8 required that values of 

C^ , Cj^ and n be available at trim. These values were generated from 

the F-8°aerodynamic data in [24] in the manner described below. 

First, by linear interpolation of moment data about the pitching axis, 

elevator deflections corresponding to trim (6 ) over a range of angle-of- 

attack (ot) were determined at selected Mach numbers. Next, C was 

L 

graphically estimated by determining the slope of a straight Ixne passed 
through points corresponding to C^ as a function of a, again at selected 
Mach numbers. Since a is small through most of a long range intercept 
trajectory, the approximation was biased to give greater accuracy to 
small a. This procedure also resulted in a non-zero ”angle-of- attack at 
zero lift" (a ) . Values of C. and a as functions of Mach number are 
shown in Table 3. The graphic^ estimates are displayed in Figures 11 
through 17. The original drag coef?xeient data is tabularized as a fxinction 
of elevator deflection, angle of attack and Mach number. The parameter 
was calculated by interpolating the drag coefficient data at trim 

conditions (6e-, a ) for various Mach numbers. Values obtained for 

1 o IJ 

are shown as a function of Mach in Table 3. Finally, n was calculated 
based on the parabolic drag model in (4.1). Values of (Cq-C^^ ) were 
plotted against (a-a^) on log-log paper and a "best fit" line^with a slope 
of 2.0 was passed through the points, again favoring lower values of a. 

Then, from values of (Cj^-Cj^ ) and (a-a^) at a point on the line, n was 
calculated by using (4.2) in (4.1) to eliminate L and solving for r\: 

o a 
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The plots used in the estimation of n are displayed In Figures 18 to 24 
and the resulting values of n are given in Table 3, 


TABLE 3 

SUMMARY OF F-8 TRIM AERO DATA 


MACH 

a© (rad) 


‘^Do 

n 

.18 

.0192 

3.366 

.0149 

.446 

.6 

.0105 

3.518 

.0142 

.580 

.85 

.0140 

4.09 

.0152 

.734 

.9 

.0143 

4.29 

.0166 

.807 

.98 

.0105 

4.18 

.0291 

.748 

1.1 

.0105 

4.24 

.040 

.733 

1.2 

.0105 

3.58 

.039 

.722 
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Figure 20. ^ a function of (a-a ) 
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at M = 0.85. 
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Section 5 


NUMERICAL RESULTS 

This section documents the numerical tests conducted to verify the 
analytical results described in Section 2 , and to test the control design 
described in Section 3. The aircraft simulated in this study was the F-8C, 
modelled as described in Section 4. The computer used in generating num- 
erical data was the IBM 370/168* This machine has a 32 bit word size, a 
basic machine cycle time of 8 ns, and performs 3 X 10^ floating point 
operations per second. Main storage access time averages 480 ns. Com- 
putation of the full zero-order solution required 0.07 seconds per update 
and 50198 bytes of core space. Control computation with proportional 
vertical lift substituted for the third and fourth boundary layer solutions 
required 0.03 seconds per update and 43724 bytes of core space. The dif- 
ference in time and space requirements is primarily due to the numerical 
search (2.43) in the third boundary layer. 

Five test cases were utilized for demonstration and evaluation of 
the control logic corresponding to a full zero-order solution. The initial 
conditions, intercept time and final range are summarized in Table 4. 

Cases 1 through 3 share the same sort of initial geometry, though with 
varying initial ranges to the target. In each case, the target flies at a 
fixed altitude, heading and velocity. These cases were chosen to display 
the effect of the range matching calculation, and hence involve very 
little maneuvering. In all three cases, the initial heading error is small 
and the aircraft is at the first boundary layer optimal altitude for its 
initial energy. Figure 25 illustrates the resulting ground tracks. All 
three cases display a slight degree of turning throughout the climb and 
descent portions of their trajectories. This is due to the fact that the 
outer solution optimal heading calculation (2.24) is based on the assumption 
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TABLE 4 


SUMMARY OF INITIAL AND FINAL CONDITIONS 
FOR TEST CASES 


Case 

x(m) 

y(m) h(m) 

8 (rad) 

V(m/s) 


1 

0. 

0. 3048. 

0. 

295.92 


2 

0. 

0. 3048. 

0. 

295.92 


3 

0. 

0. 3048. 

0. 

295.92 


4 

0. 

0. 3048. 

0. 

295.92 


5 

0. 

0. 9144. 

0. 

212.31 


Case 

Xj(m) 


Cm) 

6^ (rad) 

V^(m/s) 

1 

140208. 

-80772 . 

3048. 

1.0472 

232.56 

2 

93023. 

-53209. 

3048. 

1.0472 

232.56 

3 

46698. 

-27000. 

3048. 

1.0472 

232.56 

4 

0. 

1335.7 

3048. 

3.1416 

274.32 

5 

0. 

1335.7 

6096. 

3.1416 

274.32 


Case 

tf(s) 

r(tp(m) 

1 

618.05 

34.901 

2 

434.72 

34.008 

3 

239.98 

34.287 

4 

120.09 

17.247 

5 

150.69 

126.25 
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Figure 25. Ground tracks for cases 1-3. 




that the aircraft flies at its cruise velocity throughout its trajectory, 
ignoring the fact that lower velocities are encountered during climb and 
descent. During climb, this variation in heading is insignificant; however, 
it is more noticeable during descent. This is due to the more rapidly 
changing line-of-sight angle (X in Fig. 1) encountered as the target moves 
across the aircraft’s path at relatively close range. 

Figure 26 displays histories of actual and reference altitudes (h^) 

for Cases 1 through 3. Case 1 is a full long range intercept; the aircraft 
* * 

climbs to (E^, h^) before descending to the target. Cases 2 and 3 are 
short range intercepts in that the trajectories never achieve the optimal 
cruise energy; instead, they reach lower pseudo cruise energies prior to 
descent initiation. Note that both short range intercepts contain apparent 
cruise legs. This is because the choice of pseudo cruise energy levels was 
restricted to a set of discrete values. If a continuum of pseudo cruise 
energies was allowed, the short range trajectories would consist of 
a climb to cruise immediately followed by a descent to the target. It 
can also be seen that the aircraft altitude tends to lag behind the refer- 
ence altitude during climb in all three cases. This is because the control 
solution is based on a zero order singular perturbation analysis, which 
results in a type zero control law. Hence, a nearly constant error results 
when following the ramp- like altitude reference during climb. Inclusion of 
first-order correction terms in the control solution would be necessary to 
eliminate these errors. It is interesting to note the discontinuity at 
t = 29s in the reference altitude for Case 3. This occurred because the 
aircraft was lagging the reference altitude during climb, thus travelling 
faster than it should have at all energy levels, and closing range with 
the target at a faster pace than appropriate for the climb path. Finally, 
the range to the target became sufficiently shortened that it became 
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necessary to range match to a lower pseudo cruise energy level. A dis- 
continuity in reference altitude can also be observed at the beginning 
of descent for all three cases. This can be attributed to the energy 
state model of the aircraft and the fact that, while the aircraft’s 
optimal cruise point is defined by (2.20), the descent path is defined, 
via (2.31), as the altitude and energy corresponding to the 
In the F-8 aircraft, these points do not coincide at the pseudo crxiise 
energies encountered in Cases 1 through 3. 

Figure 27 displays aircraft altitude as a function of velocity for 
Case 2. It can be seen that the energy gain during climb occurred at 
almost constant velocity. Figure 28 displays actual and reference 
flight path angle (y^) histories for Case 2. Figure 29 shows lift and 
bank angle profiles for the same case. The beginning of descent is 
clearly marked in Figure 28 by a sudden decrease in the reference flight 
path angle. The resulting control profiles call for near maximum lift 
and inverted flight to rotate the flight path angle. This amounts to a 
zoom dive to the descent path. Figure 30 displays the thrust history for 
Case 2. During descent, ripples can be seen in both lift and thrust. 
These arise from the proportional vertical lift calculation used during 
the descent leg. The reference flight path angle, from (3.13), is 
dependent on the term E(dh^/dE)/V. The derivative is calculated by a 
first order difference expression directly from the tabular descent data 
listed in Table 2. The ripples in y^^, and hence in lift (3.15), occur 
when the aircraft altitude passes across intervals between the discrete 
points comprising the descent data. These ripples affect the thrust 
through 5T^ in (3.21) and (3.22). The average variation is a consequence 
of the range matching taking place during descent. The accuracy of this 


58 






GMD (RflD) 


O 

o 

CO 



6i 










approach is demonstrated by the small values for r(t^) in Table 4. 

Figure 31 displays the ground track for Case 4. In this case, 
the aircraft is initially at the first boundary layer optimal altitude 
for its initial energy, offset laterally from the target by a distance 
corresponding to the diameter of a turn for the aircraft at the highest 
sustainable turn rate for the initial conditions. The target moves at 
a constant altitude and velocity at a heading opposite from the controlled 
aircraft's initial heading. This case illustrates a combined initial 
turn and climb behavior, followed by a descent under near tail chase 
conditions. Figure 32 gives the actual and reference altitude histories 
for this case. Note that during the initial hard turn, the aircraft 
performs a "high speed yo-yo" maneuver, moving up in altitude (above the 
climb path) to trade speed for enhanced turning performance. During this 
phase, the reference altitude from the second boundary layer (h^) deviates 
markedly from the first boundary layer optimal altitude (h^) , due to the 
initial heading error. Within thirty seconds, however, the heading error 
has been brought down to a small value, resulting in the reference altitude 
asymptotically approaching the altitude called for from the first boundary 
layer solution. Note that the sudden change in reference altitude around 
t = 10 s is consistent with the behavior of the second boundary layer 
solution as displayed in Figure 4, where a sudden jump to the comer velo- 
city altitude takes place for large heading errors at low to midrange energy 
levels. Figure 33 shows the altitude versus velocity profile for Case 4. 
Figures 34 and 35 display the reference and actual flight path angle 
profiles, and the lift and bank angle profiles for this case. 

Figure 36 is the ground track for Case 5 . This case involved both a 
large initial heading error and a 6096 m offset from the climb path. 

This case was selected to demonstrate high altitude turning behavior, and 
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Figure 33. Altitude versus velocity for case 4. 








Figure 34. Command and actual flight path angle profile 
for case 4. 
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the balance in allocation of resources for reduction of heading and altitude 
errors. The horizontal plane geometry for this case is similar to that of 
Case 4, with the target moving in the same manner, though at a hi^er alti- 
tude. Again, the aircraft is initially separated from the target by a 
distance corresponding to the diameter of a turn made at the highest sus- 
tainable turn rate for the aircraft’s initial conditions. Figure 37 gives 
the actual and reference altitude histories for this case. It can be seen 
that initially more emphasis is given to reducing heading error. In Figure 
38, altitude is shown as a function of velocity. Figures 39 and 40 are 
the reference and actual flight path angle histories, and the lift and bank 
angle profiles. 

Comparing the lift histories for cases 4 and 5, one notes that in Case 
4 the low altitude turning behavior is characterized by a liberal use of 
maximum lift (Fig. 35), which indicates that zeroing of heading error, under 
the circumstance^ is of greater importance than gaining energy. On the 
other hand, for Case 5, lower values of lift are called for (Fig. 40). This 
is in agreement with Figure 3, which shows that the control solution for 
the second boundary layer tends to inhibit the use of lift at hi^er energies. 
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SECTION 6 


CONCLUSIONS AND RECOMMENDATIONS 

It has been demonstrated that singular perturbation methods are an 
effective tool for developing computer algorithms for on-line optimal air- 
craft controls. The computation time and storage requirements appear to be 
within the capabilities of present day flight computers. This evaluation 
is based on a minimum time intercept problem. Extensions to other performance 
indices and other forms of terminal constraints should be possible. 

An essential aspect in obtaining Implementable solutions is the ability 
to order the individual state variables on separate time scales. This re- 
quires considerable insight to the d 3 mamics of optimal flight so that a 
suitable ordering can be made at the outset. The ordering selected here 
should be applicable to a variety of problem formulations. In general, not 
all state variables will have separable dynamics. However, this study has 
illustrated two techniques for overcoming this difficulty. For short range 
intercept problems, the position and energy dynamics are coupled, and the 
problem was corrected by constraining the cruise energy level. In general, 
altitude and flight path angle dynamics are highly coupled for all intercept 
conditions. In this case, penalizing flight path angle variations is effective 
in accounting for this coupling. 

A second aspect that may be a stiimbling block in singular perturbation 
analysis (at least from the perspective of real time control) is the need to 
define the terminal boundary layer initiation times. The number of these 
layers is dependent on the number of boundary layer state variables that are 
constrained at the terminal time- This definition requires a boundary layer 
integration, which in the case of minimum time intercept can be performed 
off-line. In other formulations, boundary layer integration can have a large 
impact on the requirements for real time implementation. When they cannot 
be done off-line, then an attempt should be made to expand the necessary 
conditions to second order and obtain an analytic expression for the inte- 
gral. 

Based on the results of this study, we offer the following recommenda- 
tions : 
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1. Flight Testing - As mentioned in Section 1, the control algorithm 
has been tested at the NASA Langley real time simulation facility. 
Hence, a logical next step would be an actual flight test. Such an 
effort would require very little in terms of modifications to the 
control algorithm, and could make use of the F-8 fly-by-wire test 
bed at the NASA Dryden facility. It should be noted that, as a 
result of the piloted simulation effort at Langley, many of the 
issues relating to interfacing with displays have been addressed. 

2. Effect of Wind - To date, the effect of wind on the performance of 
the control solution has neither been evaluated nor compensated for. 
Obviously, this topic would have to be addressed before a practical 
implementation of the controller could take place. 

3. Higher Order Solutions - The control law is currently a zero-order 
approximation to the optimal control. In climb, this results in a 
significant lag in following the desired climb profile. This can be 
corrected by introducing first-order correction terms into the first 
boundary layer solution. Similarly, errors in the outer solution 
optimal heading, arising from the rapidly changing velocity and non- 
zero flight path angles in the descent boundary layer, would be 
reduced by correcting the cruise heading to first-order to account 
for these effects. These corrections would primarily affect the 
off-line computations. 

4. Other Performance Indices - Performance indices other than minimum 
time need to be Investigated. For example, a logical extension is 
to consider a weighted combination of time and fuel consumption. 

This would encompass most mission objectives. It should be noted in 
this context that the strict minimum time case is highly specialized 
in that the outer and first boundary solutions for h and V are inde- 
pendent of geometry and target velocity. 

5. Bank Angle Computation - In situations where the lift magnitude is 
small, large fluctuations result from small variations in the bank 
angle command in the horizontal and vertical lift components. In 
the context of a practical implementation, this amounts to an excess- 
ive pilot effort in return for a minor trajectory correction. Hence, 
development of a suboptimal strategy to correct this shortcoming is 
called for. 
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APPENDIX A 


OPTIMIZATION OF POSITION, ENERGY AND HEADING DYNAMICS 

This appendix provides a detailed derivation of the zero-order 
control solutions optimizing position, energy and heading dynamics * 
The state equations, for the coordinate frame shown in Figure 1, are 

(A.l) 

(A.2) 


V cos Y cos S 


y * V cos Y cos 3 - cos Yrp 


eE - (T-D)V/W 
2 


e 3* L sin y/mV cos y 
3 / 


e h* V sin y 

£ Y * (L cos u *“ WcosY)/mV 

The following conditions must be satisfied for optimality: 


(A. 3) 
(A. 4) 
(A. 5) 
(A- 6) 


= -8H/3x , Xy = -3H/3y 

(A. 7) 

eXg « -3H/3E 

(A.8) 

e^Xp- -3H/36 

(A. 9) 

-3H/3h 

(A. 10) 

e^X^- -3H/3y 

(A. 11) 

L, M, T = arg min H(x, _X, _u) 
L , y ,T 

(A. 12) 

T . 2 

H = X _x + 1 + ksinj:Y + constraints 

(A. 13) 


where the minimization in (A. 12) is subject to the constraints (2.13-2.16), 

The Hamiltonian in (A. 13) is defined for the perfomnance index in (2.12). 

A.l Outer Solution 

Taking the limit in (A.l - A. 11) as e->0 we have the following zero-order 
necessary conditions for the outer solution: 

(A. 14) 


D,y “0, Y =0, L 
o o * *o * o 


8H /3E » 3H /3h - 3H /33 * 0 
o o o 


(A. 15) 


H = X Vcos3 + X (Vsin3 - V_ cos y-) + 1 + constraints * 0 (A. 16) 
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where is drag for *» W and is defined in (2.22). Note from (A. 15) 
that E, h and 3 take the role of control variables as we let e-H) and 
(A.8 - A.ll). This implies that E, h and 3 minimize At the same 

time, the original controls now become constraints in (A.14). The optimal 
heading 3^ is derived as follows. From (A*15) we have 

3H^/33 » Vsin3 + X Vcos3 » 0 (A. 17) 

° ^o ^o 

or that 

tan3^ = X /X^ (A. 18) 

^o o 

implying, from (A. 7), that the optimal heading is constant. Referring to 
Figure 1, the geometrical requirement for intercept in the horizontal 
plane is 


Vsin(3-X) = cos cos X (A. 19) 

where X is the line-of-sight angle. 

The optimal cruise point (h^,E^) is derived in the following manner. 
We use (A. 18) to eliminate X^ from (A. 16) resulting in 

= X^ V/cos3 - X^ V^cos Y^tan 3 + 1 + constraints = 0 (A. 20) 

o ^o 

This gives 

X^ /cos3 = -l/(V-V^cosy^sin3) < 0 (A. 21) 

o 


From the intercept condition (A. 19), it can be seen that the denominator 
of (A. 21) must be positive. Applying this condition to (A. 20) and minimizing 
gives: 


h 


o* 


E 

o 


arg max(V-V^ cos y^ sin3) 
h,E 


T =D 
d o 


(A.22) 


Despite the seeming dependency of the solution of (A.22) on V^, y^ 
and 3^ it can be shown that and E^ are independent of target parameters 
and intercept geometry. Denoting the quantity being maximized in (A.22) as 
n, we will show that 3n/3V > 0 for all values of V. It follows that (A.22) 
reduces to finding h and E that maximizes V subject to the constraints (A.14) 
and also (2.13-2.16). Examining the variation of (A.22) and the intercept 
condition (A. 19) we have 
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(A.23) 


6n ■ 6V-V^cosY^cosg^ 6& 


6Vsin(0^-X) + cos (B^-X)6e 


respectively. Eliminating £3 In (A.23) and (A. 24) we have 


(A.24) 


6n/6V - 1 + V^cosy^ COS^ tan(3^-X)/V^ 


(A. 25) 


Making use of (A. 19) to eliminate in (A.25), we obtain 

6n/5V » 1 + cos3sln^(3 -X)/cos (6 -X)cos X 
o o o 


(A. 26) 


For intercept to occur, (g -X) must lie in the interval (-ir/2, r/2), implying 

o 

that cos (3^-X) > 0. Further, from Figure 1, cos3^/cosA > 0. Note also 

that if X =» 3 =• +7 t/ 2, then the second term in (A. 26) is zero, 

o — 

After determining h^, that maximize V subject to the constraint 
T^=D^, 3^ is computed from (A. 19) for the maximum cruise velocity, V=V^ 
then the position costate can be determined from (A. 18) and (A. 21): 


-cos3q/(V-V^ cos Y.p sin 3^) 


(A.27) 


-sin3^/(V-V cos si^i 2„) 
olio 


(A. 28) 


Note that (A. 18) should not be used directly in calculating X^ given X^ , 

due to the indeterminate form that results at 3 =* +7 t/2. ° ° 

o — 

A, 2 First Boundary Layer Solution 

Using the time transformation t =* t/e and letting e-K), the necessary 
conditions for the first boundary layer become: 


0, Yj^ = 0, 


(A. 29) 


3H^/3h =• 3Hj^/33 = 3H^/3T 


(A. 30) 


- X^ Vcos3 + X (Vsin3-V^cosY^) 
o ^o 

+ X_ (T-D )V/W + 1 + constraints » 0 
o 


(A. 31) 
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Although (A. 30) and (A. 31) represent four nonlinear equations in four 
unknowns, their solution is quite straightforward. 

The first boundary layer optimal heading comes from (A, 30): 

(A. 32) 


3H-/93 » -X VsinB + X Vcos3 - 0 


giving 


- 1 , 


= tan 

Because T appears linearly in (A. 31), we have 

T, 


(A. 33) 


(h,V) 

for X <0 

(A. 34) 

max 



MnO-’" 

for X >0 

^1 

(A. 35) 


The minimization with respect to h is obtained by recognizing, as 
described in Appendix E, that 9H^/9h « 0 and = 0 is equivalent to 


h = arg min 
^ h 


(T -D )V/W 
max o 

H^(E,h) 


(A. 36) 


current 


for X_ <0, and 
^1 


, . (T , -D )V/W 

h = arg min min o 

^ h H (E,h) 

o ’ 


(A. 37) 


current 


for X_ >0. The term H (E,h) is the outer solution Hamiltonian evaluated 
El o 

at the current values of E and h. The expressions (A, 36) and (A,37) roughly 
correspond to climb and descent, respectively. The first boundary layer 
costate is obtained from (A. 31): 


X = -ra (E,h)/V(T-D ) 

£il O O 


(A. 38) 


T = T, 


As was the case with the outer solution, it is possible to obtain hi 
in a form independent of target parameters. Use of the condition H^(h^,E^)*0 
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provides 


1 - X cos - -X^ 

■^O o 

Substituting (A. 39) in (A,36) and (A-37), using the inequality in (A. 21), 
and then aliminating multiplicative constant terms results in: 


^ ^ (T -D )V 

h ■ arg min max o 

^ h V - V 


current 


for X„ <0, and for X- >0: 

a- Ea~ 


h- » arg min 
h 


V - V 


E = E 


current 


(A. 40) 


(A. 41) 


It should be noted that h, (E ) , where E is a cruise or pseudo cruise 

1 o' * o 

energy level, will generally differ between climb and descent; hence the 
aircraft must execute a constant energy transition arc from cruise to 
descent altitude, as discussed in Section 3. Also of interest here is 
the fact that T- may not switch from T to T . immediately upon the 
initiation of descent. This can be seen by substituting (A. 39) in (A. 38), 
giving 



-W X (V-V )/V (T-D )cos e 
X o" ^ o' o 

o 


h = h 

1 


(A. 42) 


During climb T- > D and X_ < 0; thus, * T during climb. After the 

^ 1 o ’ 1 max ® 

initiation of descent, for large values of E^, ^nd initially, 

V > V^, so we still have that X^ <0. Because of this, the optimality con- 
dition is T- = T in descent until V < V , as illustrated in Figure 8. 

1 max o ® 

For V < V^, the optimality condition is =* ^min* Note from Figure 8 that 

if E^ < E^, then V < initially on the descent path, and thrust should 

immediately switch to T . . 

min 


A. 3 Second Bomdary Layer Solution 

The heading dynamics are accounted for by introducing the time trans- 
2 

formation t * t/e and once again letting e-K) . The necessary conditions for 
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optimality are 


0 , 


22 


(A. 43) 


3H2/3h = 3H2/3T2 * 3^2/3122 * 0 


(A. 44) 


H 2 = H^(E,h,6) + Xg L22/mV - VK L22/qsW + 
constraints » 0 


(A- 45) 


where second boundary layer horizontal lift and H^(E,h,3) is 

the first boundary layer Hamiltonian evaluated at second boundary layer 
conditions* Equations (A, 44) and (A. 45) give four nonlinear equations in 
four unknowns . 

If we assume that all turning takes place near the initial time we 
have that X < 0 so that the optimal thrust during turning is 


^2 ■ 0 ' 2 -'' 2 > 


(A. 46) 


where h .2 is the second boundary layer optimal h and is 


V - /(E-h )2g (A. 47) 

E = E 

current 

Solving for the optimal horizontal lift, using (A. 44) and (A.45), we obtain 


3H2/3L22 = -2X^ KL 22 V/qsW + X^ /mV = 0 


(A. 48) 


Using (A. 48) to eliminate X in (A.45) results in 

^2 


L 22 = /-qsWH^(E,h,e)/VKXg -sign (8^^-e) 
We also note from (A. 48) that 


(A. 49) 


sign(L 22 ) = sign(X^ ) 


(A.50) 


Substitution of these results in (A.45) and use of the minimization proce- 
dure detailed in Appendix E results in the following solution for h^: 


h 2 = arg min {-p/H^(E,h,8)KV} 
h 


8=8 


current 


(A.51) 
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subject to the constraints (A. 44) and (2.13-2.16). Having computed the 
second boimdary layer controls, the costate X is determined: 


X. =2H- (E,h,6)mV/L„ 

?2 ^ h-h2 (A. 52) 
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APPENDIX B 


OPTIMIZATION OF ALTITUDE AND FLIGHT PATH ANGLE DYNAMICS 


This appendix provides a detailed derivation of the zero-order control 
solutions optimizing altitude and flight path angle dynamics. The state 
equations and necessary conditions have already been summarized (A.1-A.13). 
B.l Third Boundary Layer Solution 

The third boundary layer problem arises from applying the time trans- 

3 

formation t * t/e to (A.l-A. 12) and letting e->0. Here, the zero-order 
necessary conditions for optimality are: 


L^3 = W cos y (B.l) 

BH^/Sy = 3H^/3L2 ° 

H3 = Kj^Ch.E.e.Y) - Xg L2^V/qsW + L2g/WV cos y 

VsinY + cos^-L^) = 0 (B.3) 

where H^(h,E,P,Y) is the first boundary layer Hamiltonian evaluated at 
current values of h,E,g and y; that is 


H^(h,E,3,Y) cosS + X sing)VcosY - X V^cos Y-j. 

o ^o ^o 

+ X (T-D )V/W + 1 + ksin^Y (B.4) 

El o 

where 

= qsCj^ + KW^cos^ y/^s (B.5) 

The parameters X , X , X_ , and X^ are known from the outer, and first 
"^o ^o ^1 ^2 

and second boundary layer solutions. Satisfaction of (B.2) for yields: 
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■dYL^/^h^ - -2Xg KV/qsW + Xg g/WVcosY - 2 V 2 L 2 - 0 


(B.6) 


Note that v»0 off the lift constraint, in which case the horizontal lift 
solution is given by 


^23 ^ ^ ^3 


(B.7) 


On the constraint botmdary, - L 

ib o max 


2 _W 2 -2 


COS Y and v can be determined 


as a function of y. In either case, and y can be determined as functions 

of Y and other known parameters. Hence, we need only combine (B.7) with 

the first of (B.2) and (B.3) to compute optimal y^ and its corresponding 

costate X, . To do this, we follow the procedure described in Appendix E, 

3 

writing 


Y^ = arg max {siny/H^ (h,E,3 ,y) } “ sign (h^-h) 
where is the second boundary layer Hamiltonian 


(B.8) 


H 2 » H^-X^ L ^3 KV/qsW + X^ 1^3 g/WVcosy 


(B.9) 


For L ^3 in (B.9), (B.7) is evaluated at the search value of y in (B. 8 ). 

The maximization in (B. 8 ) is performed over the range 0 ^ y $ '^max* 
we have used the fact that sign(X^ ) = -sign(y) - sign(h-h 2 ) where h^ is 
the optimum altitude from the second boundary layer. Having determined Y 3 
and 1«23> can be computed using (B.2): 


Xu = -{H. 


;^(h,E,B,Y3) " (B.IO) 


B.2 Fourth Boundary Layer Solution 

For solutions off the lift constraint bound we write the fourth boundary 
layer Hamiltonian as 


= HgCh.E.B.Y) - K(26Lj^^ Wcosy + + 2(6L2^L23 + L2^)V/qsW 


+ Xg 6 L 2 ^ g/WVcosY + X^ 5L^^g/WV = 0 


(B.ll) 
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where ^^24 perturbations in the controls from the third 

boundary layer solutions: 


<51-14 “ ^14 

- Wcosy 

(B.12) 

^^4 “ ^24 

- L 23 

(B.13) 


Using the conditions = 0, = 0 and 3H^/36L2^ * 0 we obtain the 


following solutions 

8qs/2Xg KV^ cosy - (B.14) 

2 1/2 

- [-H 3 (h,E,B,Y)qs W/X^ KV-BL^^] -s±g^(.y^-y) (B.15) 

The final lift and bank angle equations are: 

L = {(WcosY + + (L 23 + (B.16) 

y * tan"^(L 2 ^/Lj^^) (B.17) 


If L in (B.16) exceeds L we set L = L and reformulate the necessary 

max max 

conditions in order to determine the optimum bank angle. 


Lt , = L cos u 
14 max 


(B.18) 


L«, = L sin y 
24 max 


(B.19) 


we obtain, via 9H^/3y = 0: 

tany = X /X cos y 


(B. 20 ) 


The procedure for and considerations involved in determining X are 

^4 

detailed in Section B.3, 

It can readily be shown that the fourth boundary layer solution asymptoti- 
cally approaches that of the third. From (B.7) and (B.14), we have 


6 L 


24 


X^ qsg(l/2X_ KV^cosy - 1/2X_ KV^ cos Y-j) 

P2 c»2^ j 


(B.21) 
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Since (B.21) is continuous it will asymptotically approach zero as 
Examining the form of (B.15), it can be seen that also go to zero as 

r*yy since H^Ch.E.Sjy) goes to zero, as does 
B.3 Calculation of the Fourth Boimdary Layer Costate 

The costate X is determined by a quadratic equation obtained by 


^4 

substituting (B.20) in (B.ll) : 

AX^y, + BX + C * 0 (B.22) 

where 

A » g^(L ^/W^-cos^y)/V^ (B.23) 

max 

B * 2$g cos y/V (B.24) 

C = (Xq L g/VWcosy)^-$^ (B.25) 

max 

4. = Hj^(h,E,6,y, L=L ^^^ ) + Vsiny (B.26) 


where 


% 


has been defined in 


(B.4). 


The roots of (B.22) are as follows: 


X = (-B + y42-4AC) /2A 


(B.27) 


where the term under the radical is 


B^-4AC - 4{(«L„_g/VW)^-($L^_g/VW) - 


max 


max 


+ (L g/VW)^/cosY)^} 


(B.28) 


Because we require real— valued roots for (B.22) the term (B.28) must 
remain nonnegative, or 


+ U-(L„^^/W cosY)^} (X„ g/V)^ 5 0 

max P2 


(B.29) 
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or 


^ a ^ COS^ y) (X. g/VWcosY)^ 


max 


(B.30) 


In the third boundary layer it can be seen that 

- W^cos^Y (B.31) 

2 max 

max 

Also, the third boundary layer Hamiltonian in (B.3), evaluated at the 
conditions under consideration is 

=* $ + Xq g/VW cos Y ^ 0 (B.32) 

^ ^max 

The product X^ is always negative, so (B.30) can be rewritten: 

62 2 

5 |X- L, g/VW cos y1 (B.33) 

^2 max 

Combining (B.31) and (B.33), we can see that the inequality in (B.30) is 
guaranteed. 

It can also be shown that the roots of (B.22) have opposite sign sense. 
This is required to make a selection from the two solutions for X in (B.27) 
The importance of the sign sense of X can be appreciated by noting that 


X = 3J/3y 
^4 


(B.34) 


where J is the performance index (2.12). Since Y 3 is the optimal flight 
path angle, we have from (B.34) that 


6J * X (y-Yo) S 0 
^4 ^ 


(B.35) 


This indicates that X must have the same sign as (y“Yo) • To show that 

"^4 

(B.27) always produces roots of opposite sign, note from (B.23) that A i 0 
since ^ W within the flight envelope. Also, (B.25) and (B.33) indicate 

that C ^ 0, implying that AC $ 0 in (B.28). Because of this, the term under 
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the radical In (B. 27 ) has a dominating magnitude In that expression, 
guaranteeing real roots of opposing sign. 
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APPENDIX C 


ORDERING OF ENERGY AND HEADING DYNAMICS 


This appendix justifies the ordering selected for energy and heading 
dynamics. We will examine the eigenvalues associated with a linearized 
boundary value problem in which E and 3 dynamics are both considered on the 


same time scale (t “t/e, in the context of 
necessary conditions during climb are: 

Section 2). 

In this case, the 

dE/dt - (T -D)V/W 

max 

dXg/dx 

- -3H/3E 

(C.l) 

dS/dt - L2/mV 

dXg/dt 

= -3H/3S 

(C.2) 

3H/3h - 0 

3H/3L 

2 

- 0 

(C.3) 


H - X V cos 6 + X (V sin 6 - V_ cos y„) + X„(T -D)V/W 
Xq y q -l i Cl max 

+ XgL^/mV +1-0 (C.4) 

where X and X are defined In (2.20) through (2.26), and 
X y 

D = qsCp + K(W^ + L2)/qs (C.5) 

The second of equations (C.3) gives 

= (psg/4KX^)Xg (C.6) 

Consider an expansion of the above conditions about the optimal climb solu- 
tion at the optimal heading, so that we have 

ii - h^ (E) 3H/9h - 0 L 2 » 0 

dE/dt - E^ (E) 5 (C.7) 

All of the nominal values in (C.7) are associated with the first boundary 
layer solution in Section 2, The linearized dynamics resulting from (C.l) 
through (C.3), (C.5) and (C.6), using primes to indicate d(*)/dT, are: 
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(C.8a) 


6E' - (SE’/3E)6E + (3E'/3h)6h - (3E^/3E)6E 


6X_ - -(3fe,/3E)6X„ + f(E,h)6E 

(C.8b) 

68' - (psg/XgKmVj^)6Xg 

(C.9a) 

6X^ - (X^/cos 8)Vj_ 68 

(C.9b) 

where 6h has been eliminated by using the first equation in (C, 3) . The 
function f(E,h) does not affect the eigenvalues of (C.8). Note that (C.8) 
and (C.9) are decoupled. The eigenvalues of (C.8) are 

Pg - ±(3E^/3E) 

(C.IO) 

and for (C.9), they are 


Pg = ±y(x^/cos 8)psg/XgKm 

(C.ll) 


where (X /cos 8) and X_ come from (2.25) and (2.32), respectively; 

X £i 

Further, both quantities are always negative. It can be seen from Table 
5 that, for the problem described in this report, 3 is clearly a faster 
variable than E. Note that in the case where relative position dynamics are 
ignored [12], we have that ^ that E is the faster variable. 
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TABLE 5 


COMPARISON OF EIGENVALUES ALONG CLIMB PATH FOR AN F-4 AIRCRAFT 


Energy (M) 

PEd/s) 

Pj(l/3) 


6705.6 

± 0.0200 

± 0.317 

12192.0 

± 0.0186 

+ 0.178 

18288.0 

+ 0.0023 

± 0.164 

24384.0 

± 0.0090 

± 0.152 

29870.4 

± 0.00115 

+ 0.763 


* 

* 

30480.0 

0 

OO 


For the F-4, 3E-/3E * X * 0 in cruise due to the Mach limit. 
1 E^ 
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APPENDIX D 


SEPARATION OF ALTITUDE AND FLIGHT PATH ANGLE DYNAMICS 


Optimization of h and y dynamics is a classical problem in flight 
mechanics. This problem's difficulty lies in the fact that these dynamics 
are highly coupled if viewed separately from position and energy dynamics. 

This appendix presents summary analysis of the h and y dynamics in the con- 
text of relative position dynamics and demonstrates that, while not completely 
decoupled, they can still be optimized by singular perturbation methods. A 
more detailed development may be found in [25]. 

D.l Formulation 

In order to simplify the discussion, (2.1 - 2.6) will be specialized to 
motion in a vertical plane, with 3^ “ t^/2 and “ 0. For this problem the 
outer and first boundary layer solutions are as in Section 2 and Appendix A, 
with X *0 and X = -1/Vq. Satisfaction of boundary conditions on h and y 

^o y© 

requires further boundary layer analysis of these dynamics, using L and T as 
control variables. It has been shown [7] that if h and y are chosen on the 
same time scale: 

e^h = V sin y (D.l) 

=• (L-W cos y)/mV (D.2) 


a nonlinear TPBVP results, the solution of which is not available in feedback 
form. In hope of obtaining a valid feedback solution, we are led to con- 
sider a further separation of h and y dynamics, amending (D.2) so that y 
varies on a faster time scale than h. 


3* 

e Y 


(L-W cos y)/mV 


(D.3) 


The validity of the approximation implied in (D. 3) will be examined by 
analyzing the eigenvalues of the linearized closed loop dynamics along the 
climb path. These eigenvalues will be compared to the eigenvalues for the 
linearized closed loop dynamics resulting from the control solution when h 
and Y 3 .re modeled according to (D.l) and (D.2). It is shown that a close 
match can be obtained by properly selecting kin (2.12). 
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D.2 Eigenvalue Analysis 

We first consider an expansion of the second boundary layer necessary 
conditions for the formulation in (D.l) and (D.2). After eliminating the 
control as a function of states and adjoints, the following fourth order 
linear system results: 

0 V, 0 0 


d 

dx 


where 


6h 

Y 

6\ 


y 


gh/v^p 

K, 


0 

0 

K„ 


0 

0 

-V, 


g/V,X 
1 Y 

-gP^/V^P 


6h 

Y 

^h 

5X 


(D.4) 


h^ (E) 

» 4m X„ K/ps <0 


h 

'^E 

P = p(h) 


Ph “ Sp/3h 


K- = -3 H /3h 1 - ^ 0 
I h 

^2 = ^ ^1 - \ 8/^1 ^ 0 


Vi ./ 


H, 


2g(E-h) 


. 2 


A V + (T-D )V/W + 1 + k sin y 
y E, ' 


o 

The eigenvalues of (D.4) are the roots of 


(D.5) 


s + a s + b 


(D.6) 


where 


a = g(g/V - 2p /p) - X g/X 
b = g Vi K^/X^ 0 


(D.7) 

(D.8) 


and are arranged symmetrically about the real and imaginary axes. Since it 
is always possible to suppress two of these modes by appropriate choice of 
initial conditions on X^ and 6X^, a necessary condition for the existence of 
an asymptotically stable boundary layer solution is that none of the eigen- 
values are strictly imaginary. This follows if 


or 


4b > a 

b > 0 and a < 0 


(D.9) 

(D.IO) 
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An analysis similar to the above was conducted in [7], where relative posi- 
tion dynamics were ignored. In this case, X » 0 so that, from (D.7), 

^o 

a > 0. Thus, the conclusion regarding existence of a boundary layer solution 
depends on (D.9), which must be tested for individual aircraft types. When 
relative position dynamics are included, the last term in (D.7) is dominant, 
and a < 0. Furthermore, it can easily be shown that, for the conditions in 
(D.IO), the damping ratio is greater than 0.707. 

We now consider the boundary layer dynamics for (D.l) and (D.2), but 
for the feedback control solution of Appendix B specialized to the same 
planar problem. This implies we are using the approximation in (D.3). When 
the closed loop dynamics are linearized about the climb path, we obtain the 
following second order system: 


d (5h)/dx - Y 


d y/dx 


6L/m V 


where 


6L - Y - K, 5h 
3 4 


K 


3 - V3 + 2 Xj. KWV^/q^ s + 2k) 


■'4 ■ * \ \h> 

Kj - ./-qi S W /2 K Vj' 


1, = 3^ E/3h^ 
nn 


- p / 2 


(D.ll) 


(D.12) 


The eigenvalues of the closed loop system are given by the roots of 


s + K3 s/mV^ + K^/m - 0 


(D.13) 
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D.3 Numerical Comparison of Eigenvalues 

Table D.l displays the eigenvalues of (D.6) and (D.13) for k-0, and 
those of (D.13) for k»0. 89 for five different energy levels along the climb 
path. F-8 aircraft data was used for this comparison. Note that for k»0, 
the damping ratio obtained using (D.13) is approximately half of that ob- 
tained from (D.6). Though not apparent from Table 6, the natural fre- 
quencies in all cases are equal. 

From the second of (D.12), it is apparent that k only affects the 
damping ratio and not the natural frequency of the closed loop dynamics cor- 
responding to (D.13). Thus, for a given energy level, it is possible to 
choose k so that the eigenvalue of (D.13) match those of (D.6). For 
example, in the case of the F-8 aircraft at E * 9112m, the eigenvalues are 
matched by selecting k =* 0.89. The calculation at other energy levels is 
summarized in the last column of Table 6. Note that a reasonably good 
approximation to the eigenvalues of (D.6) is obtained at all energy levels. 

A better approximation results when k is chosen so that the eigenvalues are 
matched at E = 13528m. This approach to approximating the optimization of 
h and y is closely related to a state constrained matching approach used in 
[12]; however, the form of the solution here is much more appropriate for 
real time implementation. Trajectory results for k=0 and k=0.89 are given 
in [25]. 
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TABLE 6 


COMPARISON OF EIGENVALUES FOR h, y BOUNDARY LAYER DYNAMICS 


Energy Level (m) 

Equation(D.6) 

Eigenvalues (1/s) 
Equation (D, 13) , k=0 

Equation (D.13); k=0.89 



9112 

-.156 

± 

1 

.108 

CM 

00 

o 

r 

± 

1 .171 

-.156 

+ 

i 

.108 

11320 

-.114 

+ 

1 

.078 

-.064 

+ 

1 .123 

-.123 

± 

i 

.063 

13528 

-.089 

+ 

1 

.068 

-.048 

+ 

1 .100 

-.093 

+ 

1 

.062 

15737 

-.069 

+ 

1 

.064 

-.032 

+ 

1 .089 

-.064 

+ 

i 

.070 

17945 

-.084 

+ 

1 

.072 

-.038 

+ 

1 .105 

-.067 

+ 

i 

.090 



APPENDIX E 


MINIMIZATION OF A HAMILTONIAN FUNCTION WITH ONE UNKNOWN COSTATE 


This appendix dociunents a method for minimizing a Hamiltonian function 
with one unknown costate, the procedure having been used in minimizing the 
Hamiltonian for the outer solution and first through third boundary layers 
of the analysis in this report. 

Given the Hamiltonian function 

H(x) • f(x) + A g(x) « 0 (E.l) 


The sufficient conditions for the existence of a minimum for (E^l) are 

(E.2) 
(E.3) 


3H/3X " ^®x * ^ 


3^H/3x^ = f + Xg >0 
XX ®xx 


In a free time problem where t does not explicitly appear in H, it is also 
necessary that H*0. Using (E,l), this leads to 


X = -f/g , s * 0 (E.4) 

Using (E.4) In (E.2) and (E.3), we obtain: 

3H/3X = - (f/g)gx = 0 (E.5) 

3Vax^ = ^xx ■ ^ ° 

Define the function 

L = g/f (E.7) 

Taking the first and second partials 

3L/3x » [fg^ - gf^]/f^ (E.8) 

3^L/3x^ - [fg^^ - gf^]/f^ - 2f^(3L/3x)/f (E.9) 
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Setting (E.8) equal to zero gives the same condition as in (E.2) and (E.4). 
Condition (E.6) can be rewritten 


®^xx 

^®xx 

> 0 

> 

g 

> 0 

8^xx - 

^8xx 

< 0 

» 

g 

< 0 


(E.IO) 


Using 3L/3x = 0 in (E.9), it can be seen that the following will be equivalent 
to (E.IO): 


3\/3x^ < 0 


3^/3x^ > 0 


> g > 0 

. g < 0 

From the foregoing, we Can conclude that 

max {L} , g > 0 

min {L} , g < 0 


(E.ll) 


(E.12) 


is equivalent to the conditions in (E.l - E.3). 
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APPENDIX F 


ELIMINATION OF NUMERICAL SINGULARITIES 


This appendix describes the measures taken to prevent the appearance 
of indeterminate ratios in the control calculations, the presence of which 
led to discontinuities in the controls when the aircraft closely followed 
the optimal altitude, heading and flight path angle commands. 


The expression for the first boundary layer costate (X„ ) is, from (A. 38) 

^1 




(F.l) 


E = E 

* current 

where is the outer solution Hamiltonian and the ”1" subscripts indicate 
values at first boundary layer optimum conditions. If the aircraft’s maxi- 
mum velocity cruise point lies on the zero-energy- rate contour, an indeter- 
minate ratio results at h*h , E=E . The value of actually tends toward 

0 0 

zero as the cruise point is approached, so before calculating X_ , the energy 

^1 

rate is tested for proximity to zero. If it is too close for a reliable co- 
state calculation, X_ is set to zero. 


Both the first and second boundary layer solutions specify an optimal 
altitude (h^ and ^ 2)9 with h^ h^ as the aircraft heading error approaches 
zero. Generally, the second boundary layer optimal altitude (h^) is deter- 
mined by solving the following equation, from (A. 51), using a step search 
in altitude: 


= arg min {- p /H^(E,h, $) VK} 




current 


where 



/-q s WH^(E,h,3)/V K X^^ 


(F.2) 


(F.3) 


Since H^(E,h^,B) 0 as B (F.2) approaches an 
the control logic, aircraft heading error is tested 
If it falls within a certain region h^ is set equal 
altitude search. 


indeterminate form. In 
for proximity to zero, 
to h^, eliminating the 
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In the third boundary layer solution, an optimum flight path angle is 
calculated by solving, from (B. 8 ): 

Y 3 » arg max {sin y /H^ChjEjB.y) > * signCh^-h) (F.4) 

Y 

This is done by a step wise search in the same manner as in equation (F.2) 
Again, an indeterminate form appears in (F.4) when h h 2 » since both y^ 
and H 2 0. When [h 2 -h| is small, the current altitude is artifically 
determined to be 


h^ - h 2 + 6 h, 


6 h > 0 


(F.5) 


Requiring 5h to be positive precludes violation of the mach constraint 
boundary. Then Y 3 and are calculated in the usual manner, but for 

h=h°. Finally, since both Y 3 and X^ are zero for h=h 2 , Y 3 and X^ are cor- 
rected using a linear interpolation: 


Y^ = Y3 (h-h2)/6h 


^h^ " ^h3 


(F. 6 ) 

(F.7) 


Although all calculations in the fourth boundary layer are analytic, it 
has been found that, due to accumulated numerical inaccuracies in the preced- 
ing boundary layer solutions, there is a lack of numerical definition when 
the flight path angle error becomes small. When this occurs, the flight path 
angle is perturbed: 

Y^ “ Y 3 + sign {y- Y 3 } - 5 Y , 6 y > 0 (F. 8 ) 

The optimal horizontal and vertical incremental lifts ^^24^ 

this boundary layer are then calculated using y according to the procedure 


detailed in 

Appendix B. Then the following corrections are applied: 


«h4 - 

• 

1 (y-Y3)/'5y| 

(F.9) 


^^4 • 

1 (y-Y 3 )/<syI 

(F.IO) 
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